package com.snap.vseries.analyze;

import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.NullArgumentException;
import org.apache.commons.math3.util.FastMath;

/**
 * 2020.4.11 按扩增倍率计算后 ，
 * a为倍率基数，作为函数中的常数项，一阶导数 二阶导数和a这个常数项无关
 * b是坡度（基线斜率）,去坡度 即b=0,在此处就是无视b;
 * c是afi,按扩增倍率,c/a就等效为新的afi;
 */
public class CurveFunction implements UnivariateFunction {

    private final double c;
    private final double d;
    private final double e;
    private final double f;
    private final double g;

    CurveFunction(double... params) {
        if (params == null) {
            throw new NullArgumentException();
        }
        if (params.length != 7) {
            throw new DimensionMismatchException(params.length, 7);
        }
        this.c = params[2] / params[0];
        this.d = params[3];
        this.e = params[4];
        this.f = params[5];
        this.g = params[6];
    }

    @Override
    public double value(double x) {
        double exp1 = FastMath.exp(-d * (x - e));
        double exp2 = FastMath.exp(-f * (x - g));
        double onePlusExp1 = 1 + exp1;
        double onePlusExp2 = 1 + exp2;

        double exp1Multi2 = FastMath.exp(-2 * d * (x - e));
        double exp2Multi2 = FastMath.exp(-2 * f * (x - g));
        double awfulExp = FastMath.exp(-d * (x - e) - f * (x - g));

        double b1 = c * ((2 * d * d * exp1Multi2) / (onePlusExp1 * onePlusExp1 * onePlusExp1) - (d * d * exp1) / (onePlusExp1 * onePlusExp1)) / onePlusExp2;
        double b2 = 2 * c * d * f * awfulExp / (onePlusExp1 * onePlusExp1 * onePlusExp2 * onePlusExp2);
        double b3 = c * ((2 * f * f * exp2Multi2) / (onePlusExp2 * onePlusExp2 * onePlusExp2) - (f * f * exp2) / (onePlusExp2 * onePlusExp2)) / onePlusExp1;
        //二阶导数
        double bb = b1 + b2 + b3;
        //一阶导数
        double aa = (c * f * exp2) / (onePlusExp1 * onePlusExp2 * onePlusExp2) + (c * d * exp1) / (onePlusExp1 * onePlusExp1 * onePlusExp2);
        //返回值：曲率
        return bb / (FastMath.pow(1 + aa * aa, 1.5));
    }
}